Super-stable tomography of any linear optical device 
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Linear optical circuits of growing complexity are playing an increasing role in emerging photonic 
quantum technologies. Individual photonic devices are typically described by a unitary matrix 
containing amplitude and phase information, the characterisation of which is a key task. We present 
a constructive scheme to retrieve the unitary matrix describing an arbitrary linear optical device 
using data obtained from one-photon and two-photon ensembles. The scheme is stable on the 
arbitrarily increasable length scale of the photon packet and independent of photon loss at input 
and output ports of the device. We find a one-to-one correspondence between ideal data and unitary 
matrix, and identify the class of non-unitary matrices capable of reproducing the data. The method 
is extended for coherent state probes, which can simulate two-photon statistics with a reduced 
visibility. We analyse the performance of reconstruction to simulated noise. 



Photonic circuitry is set to play an important role in 
emerging quantum technologies Jll, connecting or imple- 
menting quantum logic gates [2j, |3j, enabling photonic 
simulators @L [B| , and coding states for quantum commu- 
nication [1,13|j with integrated optics Q enabling large 
reconfigurable circuits [9Ull|. In these contexts, each 



lossless linear optical network is ideally described by a 
unitary matrix [12| , the determination of which is a key 
task. While a general quantum optical process may be 
fully characterised with coherent light and homodyne de- 
tection [l]| , the device to be determined and the probing 
set-up must be stable with respect to one another on a 
sub- wavelength scale. 

Since the phases imparted to each photon of a two- 
photon state combine to a single glob al p hase, circuit 
probing via two-photon interferometry |14l4l7| is stable 
on the scale of the photon packet, which is typically or- 
ders of magnitude greater than wavelength. Photons 
generated from spontaneous parametric down conversion 
(SPDC) can be lengthened through spectral filtering or 
with a cavity in narrowband SPDC [la ], while photons 
generated in an atom-cavity system can be many meters 
long fl9| . Two-photon sources are commonplace in quan- 
tum optics labs, with quantum interference visibilities 
approaching unity (20| . and future many-photon experi- 
ments are likely to be easily configured to produce sub- 
sets of data from one-photon and two-photon ensembles 
for characterisation purposes. Additionally, two-photon 
quantum statistics can be classically simulated by two co- 
herent beams with a mutually randomised phase pi. 22| 



so that any circuit information retrievable with a two- 
photon state is also available to classical states. 

In Ref. [lfl the transfer matrix of a four-mode linear 
optical device was numerically approximated as the ma- 
trix that minimised the distance between a self-simulated 
data-set of all possible two-photon statistics and experi- 
mental data, with data from one-photon ensembles used 
as starting parameters. However, a numerical search 
among the space ofmxm dimensional unitary matrices 
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FIG. 1: Circuit tomography with one-photon and two- 
photon probes is insensitive to phase fluctuations and inde- 
pendent of loss at input and output ports. 

siblc two-photon statistics becomes challenging for large 
systems. More fundamentally, it has not been clear up to 
this point if data from one-photon ensembles (one-photon 
data) and data two-photon ensembles (two-photon data) 
are sufficient to uniquely determine the matrix descrip- 
tion of a linear optical circuit. 

Here, we give a constructive method to efficiently de- 
termine the set of matrices that produce a given set of 
ideal one-photon and two-photon data, for an arbitrary 
linear optical network. Reconstruction is insensitive to 
fluctuating phases between input ports and between out- 
put ports, as shown in Fig [TJ and therefore does not 
require active stabilisation or other techniques for sta- 
ble classical interferometry. In addition, reconstruction 
is independent of the typically unknown and tomograph- 
ically troublesome photon losses at device input and out- 
put ports. The scheme involves injecting one-photon and 
two-photon states into the linear optical network with 
correlated photon detection, as shown in Fig [TJ Matrix 
amplitude information is probed with one-photon ensem- 
bles, while specific phase parameters are targeted by spe- 
cific two-photon ensembles, so that the required number 
of measurements scales quadratically with the number of 
modes, and linearly with the number of unknown param- 
eters. Assuming the device can be described by a unitary 
matrix, applying 2m of the m 2 orthonormal constraints 
leads to a one-to-one mapping between idealised exper- 
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imental data and unitary matrix, up to the fluctuating 
phases. If the port-losses are known to be zero, then a 
matrix can be recovered with no orthonormal constraints 
invoked or required. We consider non-ideal situations 
and simulate noisy data from circuits of up to 20 optical 
modes to quantify the fidelity of reconstruction versus 
noise. 

An insightful example of the tools used in our method 
is the case of a two-mode device, or subsection of a larger 
device, as shown in Fig[2ja). Sk photons are individually 
injected into input port k, coupled with efficiency 
transmitted to output port j with probability Tj k , and 
detected with efficiency Tj, giving a count rate of Rj. k 
27 |. With Rj.k = Sk Sk Tj t k fj, four such measurements 
combine to constrain four transition probabilities {Tj^} 
independently of coupling efficiencies: 
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In the case that the 2x2 device behaves unitarily, prob- 
ability constraints (e.g. T h i = T 2 , 2 = 1 - T 1)2 = 1 - T 2 ,i) 
allow description in terms of one parameter, 



rp2 
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to derive the experimentalist's loss- insensitive equation 
for the reflectivity of a beamsplitter, T^x = Br, 



Bu = 



where X replaces the RHS of (Q}. 

An alternative way to find the reflectivity of a beam- 
splitter uses the visibility of quantum versus classical cor- 
relations in a Hong Ou Mandel (HOM) dip Q. More 
generally, the visibility of a HOM dip can be used to ap- 
ply loss-insensitive constraints to the 2x2 submatrix S2 
that describes a two mode system, which may be part 
of a larger device. S2 consists of real positive probabil- 
ity amplitudes r and unit complex phase factors e za ; pre 
and post multiplication by diagonal matrices D c and D ( i 
includes coupling and detection efficiencies, to describe 
of the overall subprocess P 2 = Dd S2 D c , 
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as shown in Fig[2Jb). 

Simultaneously injecting two photons into the P 2 de- 
vice and counting subsequent coincidental detector clicks 
as a proportion, gives a quantum signal equal to the 



1.U.2} 



2,{1,2} 




■/Ws h | 



FIG. 2: Probing of circuit sub-section, (a) One-photon prob- 
ing of a two- mode subsection of larger device, (b) Two-photon 
targeting of phase parameters, with R c indicating correlated 
detection. 

square of the absolute value of the permanent of P 2 : 

Qg,hj,k = \Per(P 2 )\ 2 

12 2,2 2 \ 

— Sh Sk T g Tj [Tj k T g h -(- T g k Tj h ) 
+ S h S k r g Tj Tj ik Tj th T g ^ k Tg t h 

x 2cos(a iife - otj t h - a g ,k + a g . h ). (2) 

Introducing a delay between the photons so that no quan- 
tum interference occurs gives the classical signal equal to 
the square of the permanent of the matrix containing 
absolute values of the entries in P 2 : 
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(3) 



and can be calculated directly from one-photon data. 
The loss insensitive visibility is calculated as 



g,h,j,k 
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C g> h,j,k 
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(4) 
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Probability amplitudes r are the square root of the 
transition probabilities T described in ([IJ so can be con- 
strained from 
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A loss-independent measurable constraint on S2 phase 
angles then follows from 



Ug,h,j,k — %g,h,j,k + x g ,h,j,k 

COS (Oij.k - Oijfr — Ug.k + (Xg h) 



Vg t h,j,k Vg.hJ.k 



(6) 



Since cosine is a symmetric function, the sign of the ar- 
gument is unknown. However, the sign for any one of 
the four phase angles, say a Qt h, can be ascertained from 
([6]) if its absolute value < \a g ,h\ < T is known, and the 
absolute values and signs of the other three (reference) 
angles are known, which will be the case, as we will see. 
The sign on a gj h can be found by labelling the argument 
from © as P g ,h,j,k = \otj,k ~ a jM - a g , k + a g ,h\ then 
computing 

sgn[a g ,ft] = sgn[ 

\Pg,h,j,k - \0ij,k - Otj,h - <Xg,k - \otg,h\\\ 
-\fig,h,j,k — \otj,k — <Xj,h - Otg,k + Wg,h\\\ 1 ( 7 ) 

assuming the 3 reference angles do not sum to zero [Hj]. 

Using the amplitude and phase constraints of (J5J) and 
(|6]) with the sign disambiguation of ((TJ), we can recon- 
struct M, the m dimensional matrix description of an 
arbitrary m mode optical network, up to 2m — 1 of the 
sa m 2 parameters, independently of any unknown losses 
on input and output ports. When M can be assumed 
to behave as a unitary transformation, constraints give 
a tractable solution system for the 2m — 1 unknown pa- 
rameters. 

We can specify some equivalence among matrices 
where corresponding linear optical networks within the 
equivalence class produce the same data. Firstly, two 
matrices M a and Mb are equivalent if there exist two 
diagonal unitary matrices and such that M a = 
DY MbD^ ■ In a linear optical network M, these diagonal 
matrices can be regarded as the typically unknown and 
somewhat trivial phases on the input and output ports of 
the network, to which our data is insensitive, as long as 
photons are input in a Fock state (not superposed across 
modes) and there is no post-M interferometry before de- 
tection. These trivial phases may equally be regarded as 
a particular instance of the fluctuating phases from Fig [I] 
We encapsulate this idea in the standard way by insist- 
ing that our M is real bordered, that is, the entries on the 
left-most column and uppermost row should be real and 
positive: {Mj t i, Mi,fe}™fc = i G 1R + - Secondly, noting the 
equality of photon statistics under the complex conjugate 
transformation M — > M*, we insist that the imaginary 
part of the element be non negative: Im(M2, 2 ) > 0. 

Additionally, we make the following assumptions: for 
a straightforward description, we firstly assume that the 
two leftmost columns and two uppermost rows contain 
no vanishingly small elements — if this were not the 



case, one-photon data would identify zero matrix ele- 
ments and the reconstruction method can be straight- 
forwardly adapted; secondly, while our method is insen- 
sitive to loss on the input and output ports, we assume 
any losses are not near-total. With these equivalencies 
and assumptions, the matrix to be recovered can be pa- 
rameterised as 



M 1 



Tl,2 



n, 



t 2 ,i r 2 , 2 e' 



\T m ,l Tm^e 1 



(8) 
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where, as before, the r correspond to transition ampli- 
tudes and the a correspond to phase angles. 

Replacements are made on all elements of Mi apart 
from the first column and first row. The r amplitudes are 
replaced following ([5]) fixing Tj^ = while the indices 
of T 3j h are taken from the r to be replaced. All r not 
in the first column or row, become functions of the r in 
the first row and column. Phase angles — tt < a < tt are 
found (up to a sign) following ^ with the same index 
pattern used for the r, so that only one angle appears in 
the cosine. 

Apart from the sign on a2,2, which is defined to be 
positive, all signs on the a in Mi are determined using 
(|TJ) . Signs on the a in the 2nd column are found by fixing 
otj t k = 0-2,1 while the indices of a g< h — {a3... m ,2}, moving 
from third row to row m. Similarly, signs on the 2nd row 
arc found by fixing a jt k = ai,2 and a g . h = {a 2 ,3...m}j 
moving from third column to column m. Signs on the 
remaining a are found by fixing aj & = 0^2 and the in- 
dices of oig^h are taken from the a to be replaced. This 
pattern ensures that signs on all references phases in ([7]) 
are known, as is required. With these replacements, the 
transfer matrix can be parameterised by 2m — 1 unknown 
border amplitudes 



A/2 



/ r M T l,2 
72,1 £2,2,1,16' 

V r m,l £ TO ,2,l,ie 



ia m , 2 



••• ^2,m,l,l e " 
%m.m,l , 1 6 



where, for ease of notation, we have introduced 

Tj,hTg,k 



%g,h,j,k — %g 7 h,j,k 



T j,k 



1 



a g . h = sgnfag.h] arccos(- V g ^,j,k y g ,h,j,k), 

highlighting that the a are fully determined from exper- 
imental measurements and not functions of the r. 
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The matrix E, describing the complete experimental 
set up comprising the device M and losses on input and 
output ports, has transition probabilities corresponding 
to the set of one photon count rates R, while phase angles 
are the calculated a: E is therefore fully characterised 
so that the output state can be calculated for any ar- 
bitrary photon number Fock state input. This implies 
that probing the system with more than two photons, 
or using number resolving detectors, provides no extra 
information to solve the remaining 2m — 1 unknown pa- 
rameters in M%. In the case of no losses M = E and full 
reconstruction has been achieved. However, zero losses 
on input and output ports is atypical, particularly when 
one considers collection efficiencies from a photon source 
and detection efficiencies. 

When M can be assumed to behave unitarily, normal- 
isation can be applied to the first column of M2 and or- 
thogonality constraints applied between the first column 
and the remaining m— 1 columns, with the same applying 
for rows. The solution system is the matrix of coefficients 
in M 2 , 



6 4 (% noise for m=4) 



n 
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applied to the transition probabilities, to give the or- 
thonormal constraints, 



Aft 



w 



'1,2 



(9) 



A standard method for solving systems of linear equa- 
tions may then be used [29j |. 

One-photon and two-photon statistics can be esti- 
mated using classical coherent states. The classical ana- 
logue of one-photon counts Rj t k are intensities Ij t k, mea- 
sured at port j from injecting a coherent beam of inten- 
sity Iq into port k: the one- photon counts in (TTJ) can be 
replaced by a ratio of classical intensities. The classical 
counterpart of Hong Ou Mandel correlations are those of 
Hanbury-Brown and Twiss [231 ] arising from interference 



between two independent fields [24j, |25[ , which have al- 
ready been exploited in linear optical devices to simulate 
two- photon statistics [U [22j ■ The quantum coincidence 
rate Q g ,hj,k can be accurately estimated by injecting two 
coherent states with the same mean photon number and 
a randomised relative phase into respective input ports h 




0.05 0.1 0.15 0.2 
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FIG. 3: Average fidelities of reconstructed unitary matrices 
versus simulated noise for m=4 (top curve) and m=20 (bot- 
tom curve). 

and k. The intensity correlation function, T g ^,j,k, mea- 
sured at output ports g and j is 

T gMk = II \Per (P 2 )| 2 + / + /,./. /,,.-, (10) 

so that the quantum correlation can be estimated from 
r ' g ,h,j,k minus the product of individual intensities 

I Qg,h,j,k = Fg,h,j,k - Ij.klg.k ~ Ij,h.Ig,h- (H) 

Experimental noise and device imperfections will drive 
the reconstructed matrix away from that which best de- 
scribes the optical device, and away from unitarity since 
we use only some of those constraints when solving for 
border amplitudes. The unitary matrix U p closest to the 
reconstructed matrix may be found from the polar de- 
composition (2(|, but in the presence of errors U p will 
not, in general, perfectly match any true underlying uni- 
tary matrix Ut- 

We analysed the effect of noisy data on our reconstruc- 
tion procedure, for devices from dimension m = 4 to 
m = 20, over a range of noise level S. Noise was added to 
the simulated experimental values for one-photon count 
rates R and quantum interference visibilities V: the ideal 
values were multiplied by (1 + e), where e was chosen from 
a normal distribution centred at 0. The amount of noise 
S was controlled by varying the width of the distribution 
with S = 3<t, the third standard deviation, making it very 
unlikely that experimental values are found to be much 
further than S away from their ideal values. For a partic- 
ular value of m and <5, we generated 1000 Haar random 
unitary matrices, each of which was pre and post mul- 
tiplied by diagonal matrices with real random values to 
simulate loss. The algorithm was ran for each instance, 
with a unitary matrix extracted via the polar decompo- 
sition. The average trace distance was computed over 
all 1000 examples to give an average fidelity of recon- 
struction; the variance of each 1000 example set was also 
computed. 

Over the tested range of matrix size m, we found the 
empirical relation between average fidelity of reconstruc- 
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tion F and error level S 



F w cxp(-A(5 2 ), 



(12) 



with A m ~ 3 giving good predictions for fidelities down 
to at least sa 85%. Figure |3] shows 51 simulated data 
points between 5 = 0% and S = 5% for m = 4, and 
between S = 0% and S = 0.25% for m = 20, with error 
bars on each point corresponding to the variance of its 
1000 example data set. The empirical curve F is shown 
as the solid line. 

The effect of errors analysed here involved no attempt 
at error correction. A useful line of research would be to 
investigate how the impact of errors may be minimised 
by advantageously permuting the matrix prior to recon- 
struction. Additionally, having used only m 2m 2 two- 
photon measurements from the possible w m 4 , it would 
be interesting to understand how extra subsets of two- 
photon data may increase robustness to noise, or whether 
resources should be directed toward reducing errors on 
the minimal set. 

We have presented a scheme in which the unitary ma- 
trix of a linear optical device may be constructively re- 
trieved without the requirement for stable classical inter- 
fcrometry. The merits of probing with coherent states 
versus single photon states is interesting. Since the co- 
incidental signal falls off quadratically (due to the num- 
ber of ways of choosing 2 from m modes) an arbitrar- 
ily bright coherent state may seem preferable. However, 
two-photon states readily available within a larger multi- 
photon experiment may exhibit the same characteristics 
as those producing the larger data set, with statistics 
more representative of a device that is somewhat de- 
fined with respect to its photons; for example, the re- 
flectivity of a waveguide directional coupler is a func- 
tion of wavelength. Any advantage from an increase in 
signal and greater measurement precision with classical 
states should be weighed against potential inaccuracies 
if the crucial properties of those states do not match the 
properties of the photons, due to some experimental un- 
certainty. Furthermore, the capability to re-take tomo- 
graphic data as a subset of principal experimental data, 
for a freshly reconfigured system, during the course of an 
uninterrupted experiment, is appealing. 
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Appendix 1: A four mode example 

We present a four mode example of simulated data E < - Slm ' > for a system comprised of a device described by a Haar 
random unitary U^ 1 ™ 1 ^ with losses on inputs l/ m ) and outputs 



j^(out)jj-(sim) jr (in) 



/0.46 0. 0. 0. \ /0.245 0.54 

0. 0.65 0. 0. 0.492 0.377 + 0.192* 

0. 0. 0.41 0. 0.662 0.007 + 0.119* 

\ 0. 0. 0. 0.37/ \0.509 -0.633 - 0.34*- 



0.537 
-0.634 + 0.213*' 
0.305 - 0.339*' 
-0.042 + 0.235* 



/0.009 0.204 

0.026 0.201 + 0.102*- 
0.022 0.002 + 0.04*' 

\0.015 -0.192- 0.103*- 



0.136 
-0.227 + 0.076* 
0.069 - 0.076* 
-0.009 + 0.048*' 



0.066 \ 
0.004 -0.057*' 
-0.054 + 0.019*' 
0.035 + 0.008* / 



0.601 \ 
0.027- 0.362* 
-0.549 + 0.196*' 
0.398 + 0.095* / 



/0.08 0. 0. 0. \ 

0. 0.82 0. 0. 

0. 0. 0.55 0. 

\ 0. 0. 0. 0.24/ 



Following ^ to ([7]) with perfectly simulated quantum interference data and power ratio measurements, the unitary 
parameterised by 7 probability amplitudes is 



M. 



(sim) 



T 2,l 
r 3,l 
\ T 4,1 



(0.348 1 
(0.005 4 
-(0.565 



n,2 

0.177*')^fi 
0.082*')^^ 
_0.304*')^ki 



-(0.588 
(0.211 - 
-(0.038 



n,3 

- 0.197*)2i2ZM 
0.234*) 
-0.211*') 



Tl.l 

Tl,3T4 | i 



(0.022- 
-(0.338 
(0.319- 



n,4 

0.300*)^^- 
-0.121*)^^i 
0.076*)^^- 



giving the solution systems for the unknown r as 



/i i i 

1 0.348 + 0.177*' 0.005 + 0.082*' 

1 -0.588 + 0.197*- 0.211 -0.234? 

\1 0.022 - 0.300* -0.338 + 0.121*- 



1 \ 


/0.245 2 ^ 
0.492 2 






-0.565 - 0.304*' 







-0.038 + 0.211*- 


0.662 2 







0.319 + 0.076*- / 


\0.509 2 / 




w 





i 




1 






1 






fQ.2A5 2 \ 
0.540 2 




('\ 


1 


0.348 - 


0.177*' 


-0.588- 


-0.197*' 


0.022- 


f 0.300* 











1 


0.005- 


0.082*' 


0.211 + 


0.234*- 


-0.338 


-0.121 


i 




0.537 2 







V 


-0.5651 


- 0.304*- 


-0.038- 


- 0.211*- 


0.319- 


- 0.076* 


J 




^0.601 2 / 




w 



in agreement with £/( s ™) . 



